Numerical Analysis 



You could say that some of the equations that you encounter in describing physical systems can't be 
solved in terms of familiar functions and that they require numerical calculations to solve. It would be 
misleading to say this however, because the reality is quite the opposite. Most of the equations that 
describe the real world are sufficiently complex that your only hope of solving them is to use numerical 
methods. The simple equations that you find in introductory texts are there because they can be solved 
in terms of elementary calculations. When you start to add reality, you quickly reach a point at which 
no amount of clever analytical ability will get you a solution. That becomes the subject of this chapter. 
In all of the examples that I present I'm just showing you a taste of the subject, but I hope that you 
will see the essential ideas of how to extend the concepts. 

11.1 Interpolation 

Given equally spaced tabulated data, the problem is to find a value between the tabulated points, and 
to estimate the error in doing so. As a first example, to find a value midway between given points use 
a linear interpolation: 

f(x + h/2)^±[f(x ) + f(x + h)] 

This gives no hint of the error. To compute an error estimate, it is convenient to transform the variables 
so that this equation reads 

f(v)*\[f(k)+n-k)] : 

where the interval between data points is now 2k. Use a power series expansion of / to find an estimate 
of the error. 

f(k) = f(0) + kf'(0)+ 1 -k 2 f"(0) + ..- 
f(-k) = /(0) - kf'(0) + ^k 2 f"(0) H 

Then 

^ [/(*) + /(-*)] «/(o)+[^7"(o)], (11.1) 

where the last term is your error estimate: 

Error = Estimate - Exact = +k 2 f"{0)/2 = +h 2 f '(0)/8 

And the relative error is (Estimate — Exact)/Exact. 

As an example, interpolate the function f{x) = 2 X between and 1. Here h = 1. 

2 V2 ~ 1 [2° + 2 1 ] =1.5 

The error term is 

error » (In 2) 2 2 x /8 for x = .5 
= (.693) 2 (1.5)/8 = .090, 

and of course the true error is 1.5 — 1.414 = .086 

You can write a more general interpolation method for an arbitrary point between Xq and Xq + Ji. 
The solution is a simple extension of the above result. 
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V-fo = (x-x )(fi-f )/h, 




X Q X! 



fo = f(x ), fi = f(x + h). 
At the point x = Xq + ph you have 

V = fo + (ph)(fi - f )/h = fo(l -p) + fiP 

As before, this approach doesn't suggest the error, but again, the Taylor series allows you to work it 
out to be [h 2 p(l -p)f"(x +ph)/2\. 

The use of only two points to do an interpolation ignores the data available in the rest of the 
table. By using more points, you can greatly improve the accuracy. The simplest example of this 
method is the 4-point interpolation to find the function halfway between the data points. Again, the 
independent variable has an increment h = 2k, so the problem can be stated as one of finding the 
value of /(0) given f(±k) and f(±3k). 
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f(k) = /(o) + kf(o) + h 2 f"(o) + ^ 3 /"'(o) + • • • 

I want to isolate /(0) from this, so take 

f(k) + f(-k) = 2/(0) + k 2 f" (0) + ^A; 4 f "(0) + • • • 
/ (3 fc) + f(-3k) = 2/(0) + 9A; 2 A0) + ^k*f"(0) + ■ 
The biggest term after the /(0) is in k 2 f"(0), so I'll eliminate this. 

[/(3fc) + f(-3k)] -9[f(k) - f(-k)] « -16/(0) + 
1 



81 _ 9 
12 ~ 12 



fc 4 r"(0) 



/(o) 



16 



[ - /(-3fc) + 9f(-k) + 9f(k) - f(3k)] - [ - ^/""(0)] . 



The error estimate is then -3/i 4 / //// (0)/l28. 

To apply this, take the same example as before, f(x) = 2 X at x = .5 



2 i/2 



16 



[-2" 1 + 9-2° + 9-2 1 -2 2 



45 
32 



1.40625, 



(11.2) 



(11.3) 



and the error is 1.40625 — 1.41421 = —.008, a tenfold improvement over the previous interpolation de- 
spite the fact that the function changes markedly in this interval and you shouldn't expect interpolation 
to work very well here. 
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11.2 Solving equations 

Example: sinx — x/2 = 

From the first graph, the equation clearly has three real solutions, but finding them is the problem. 
The first method for solving f(x) = is Newton's method. The basic idea is that over a small enough 
region, everything is more or less linear. This isn't true of course, so don't be surprised that this method 
doesn't always work. But then, nothing always works. 




A general picture of a function with a root is the second graph. In that case, observe that if x$ 
is a first approximation to the root of /, the straight line tangent to the curve can be used to calculate 
an improved approximation. The equation of this line is 

y-f(xo) = f'(xo)(x-x ) 

The root of this line is defined by y = 0, with solution 

x = x - f(x )/f(x ) 

Call this solution x±. You can use this in an iterative procedure to find 

x 2 = x 1 -f(x 1 )/f'(x 1 ), (11.4) 

and in turn X3 is defined in terms of X2 etc. 

Example: Solve sinx — x/2 = 0. From the graph, a plausible guess for a root is Xq = 2. 

Xi = Xo — (sin Xo — Xo/2)/ (cos Xq — 1/2) 

= 1.900995594 f{x 1 ) = .00452 

X2 = X\ — (sin X\ — Xi/2)/ (cos x± — 1/2) 

= 1.895511645 f(x 2 ) = -.000014 

X3 = x 2 - {s'mx 2 -x 2 /2)/(cosx 2 - 1/2) 

= 1.895494267 f(x 3 ) = 2 x 10~ 10 

Such iterative procedures are ideal for use on a computer, but use them with caution, as a simple 
example shows: 

f(x) =x 1 / 3 . 

Instead of the root x = 0, the iterations in this first graph carry the supposed solution infinitely far 
away. This happens here because the higher derivatives neglected in the straight line approximation are 
large near the root. 
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A milder form of non-convergence can occur if at the root the curvature changes sign and is 
large, as in the second graph. This can lead to a limit cycle where the iteration simply oscillates from 
one side of the root to the other without going anywhere. As I said earlier this doesn't always work. 

A non-graphical derivation of this method starts from a Taylor series: If Zq is an approximate 
root and Zq + e is a presumed exact root, then 

/(2b + e) = = f(zo) + ef(zo) + ■ ■ ■ 

Neglecting higher terms then, 

e = -f(z )/f(z ), and z\ = z + e = z - f (z ) / f (z ) , (11.5) 

as before. Here z appears instead of x to remind you that this method is just as valid for complex 
functions as for real ones (and has as many pitfalls). 

There is a simple variation on this method that can be used to speed convergence where it is 
poor or to bring about convergence where the technique would otherwise break down. 

xi = x -wf(x )/ f'(xo) (11.6) 

W is a factor that can be chosen greater than one to increase the correction or less than one to decrease 
it. Which one to do is more an art than a science (1.5 and 0.5 are common choices). You can easily 
verify that any choice of w between and 2/3 will cause convergence for the solution of rr 1 / 3 = 0. You 
can also try this method on the solution of f(x) = x 2 = 0. A straight-forward iteration will certainly 
converge, but with painful slowness. The choice of w > 1 improves this considerably. 

When Newton's method works well, it will typically double the number of significant figures at 
each iteration. 

A drawback to this method is that it requires knowledge of f'(x), and that may not be simple. 
An alternate approach that avoids this starts from the picture in which a secant through the curve is 
used in place of a tangent at a point. 

Given f(xi) and f(x 2 ), construct a straight line 



y ~ f(x 2 ) = 




This has its root at y = 0, or 

*=**-f^ f(x Z~-?( Xl ) <1L7) 

This root is taken as x% and the method is iterated, substituting x 2 and X3 for x\ and x 2 . As with 
Newton's method, when it works, it works very well, but you must look out for the same type of 
non-convergence problems. This is called the secant method. 

11.3 Differentiation 

Given tabular or experimental data, how can you compute its derivative? 

Approximating the tangent by a secant, a good estimate for the derivative of / at the midpoint 

of the (x\,x 2 ) interval is 



[f(x 2 )-f(x 1 )]/(x 2 -x 1 ) 
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As usual, the geometric approach doesn't indicate the size of the error, so it's back to Taylor's series. 
Given data at points x = 0, ±h, ±2h I want the derivative /'(0). 

f(h) = /(o) + hf (o) + \tff' (o) + h 3 f"(o) + ■■■ 

fi-h) = /(o) - /i/'(o) + \h 2 r (o) - ^ 3 r(o) + • • • 

In order to isolate the term in f'(0), it's necessary to eliminate the larger term /(0), so subtract: 

f(h) - f(-h) = 2hf'(0) + \h 3 f"'{0) + • • • , 

(11 8) 

giving ^ [f(h) - f(-h)] » /'(0) + [;U 2 /"'(0) n 

and the last term, in brackets, estimates the error in the straight line approximation. 

The most obvious point about this error term is that it varies as h 2 , and so indicates by how 
much the error should decrease as you decrease the size of the interval. (How to estimate the factor 
/"'(0) I'll come to presently.) This method evaluates the derivative at one of the data points; you 
can make it more accurate if you evaluate it between the points, so that the distance from where the 
derivative is being taken to where the data is available is smaller. As before, let h = 2k, then 

^[f(k)-f(-k)]=f'(o) + h 2 f"(o) + ..., 

or, in terms of h with a shifted origin, 

l[f(h)-f(0)} -/W2) + ^ 2 r(|), (H.9) 

and the error is only 1/4 as big. 

As with interpolation methods, you can gain accuracy by going to higher order in the Taylor 

series, 

f(h) - f(-h) = 2hf(o) + ^ 3 r'(o) + ±h 5 r (o) + ■ ■ ■ 

f(2h) - f(-2h) = 4hf'(0) + ^ 3 r'(0) + ^h 5 f v (0) + • • • 

To eliminate the largest source of error, the h 3 term, multiply the first equation by 8 and subtract the 
second. 

24 r- 

8[f(h) - f(-h)] - [f(2h) - f(-2h)] = 12/1/(0) - -h 5 f v (0) + • • • , 

or 

/'(°) « T7a [f(~ 2h ) ~ 8 f(~ h ) + 8 /W " f( 2h )] - [ - ^ h 'f V (°)] (11-10) 



12h 

with an error term of order h 4 . 

As an example of this method, let f(x) = s'mx and evaluate the derivative at x = 0.2 by the 
2-point formula and the 4-point formula with h=0.1: 

2-point: — [0.2955202 - 0.0998334] = 0.9784340 

U. £ 

4-point: ^ [0.0 - 8 x 0.0998334 + 8 x 0.2955202 - 0.3894183] 

= 0.9800633 
cos 0.2 = 0.9800666 



11 — Numerical Analysis 



Again, you have a more accurate formula by evaluating the derivative between the data points: 



h = 2k 



f(k) - f(-k) = 2*/'(0) + ^ 3 f"(0) + ±k*r (0) 

97 943 
/(3fc) - /(-3fc) = 6fc/'(0) + y fc 3 /"'(0) + ~^k 5 r (0) 

91 fi 

27[/(fc) - /(-*)] - [/(3fc) - f(-3k)] = 48k f (0) - ^A; 5 / u (0) 
Changing to h/2 and translating the origin gives 

^ [/(-/i) - 27/(0) + 27 f(h) - f(2h)] = f{h/2) - A_h 4 r(h/2), 

and the coefficient of the error term is much smaller. 

The previous example of the derivative of sin a; at x = 0.2 with h = 0.1 gives, using this formula: 

1 



(11.11) 



2.4 



[0.0499792 - 27 x 0.1494381 + 27 x 0.2474040 - 0.3428978] = 0.9800661, 



and the error is less by a factor of about 7. 

You can find higher derivatives the same way. 

f(h) = /(0) + hf'(0) + \h 2 f"(0) + ^ 3 /"'(0) + ^ 4 f "(0) 

f(h) + f(-h) = 2/(0) + h*f"(0) + ±h*f""(0) + ■ ■ ■ 

r(0) = /(-*)- 2/(0) + /(*) -^(0) + ... 

Notice that the numerical approximation for /"(0) is even in h because the second derivative is un- 
changed if a: is changed to —x. 

You can get any of these expressions for higher derivatives recursively, though finding the error 
estimates requires the series method. The above expression for /"(0) can be viewed as a combination 
of first derivative formulas: 

f"(0)n[f'(h/2)-f'(-h/2)]/h 



(11.12) 



f(h) - /(0) /(0) - f(-h) 



h h 
[f(h)-2f(0) + f(-h)]/h 2 



(11.13) 



Similarly, the third and higher derivatives can be computed. The numbers that appear in these numerical 
derivatives are simply the binomial coefficients, Eq. (2.18). 

11.4 Integration 

The basic definition of an integral is the limit of a sum as the intervals approach zero. 




J2f(&)( x i+l- x i) (Xi<ii<X i+l ), (11.14) 



6 6 £3 £4 £5 
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This is more fully explained in section 1.6, and it is the basis for the various methods of numerical 
evaluations of integrals. 

The simplest choices to evaluate the integral of f(x) over the domain Xq to Xq + h would be to 
take the position of £ at one of the endpoints or maybe in the middle (here h is small). 

xo+h 

f(x)dxfaf(x )h (a) 

x 

or f(x + h)h (b) (11.15) 

or f(x + h/2)h (midpoint rule) (c) 

or [f(xo) + f(xo + hj] h/2 (trapezoidal rule) (d) 

The last expression is the average of the first two. 

I can now compare the errors in all of these approximations. Set Xq = 0. 



j\xf{x) = j dx[f(0) + xf'(0) + \x 2 f"(0) + ^ 3 /"'(0) 



+ 



hf(0) + ^/'(O) + \h^f" (0) + ±h'f"(0) + • • 



This immediately gives the error in formula (a): 

fh 



error (a) = hf(0) - J dxf(x) « -\h 2 f'{Q). (11.16) 



The error for expression (b) requires another expansion, 

fh 



error 



(b) = hf(h)- [ dxf(x) 
Jo 



= h [/(o) + hf (o) + ■■■]- [hf(o) + 1/iV'(o) + • • • ] 

« ^ 2 /'(0) (U.17) 

Since this is the opposite sign from the previous error, it is immediately clear that the error in (d) will 
be less, because (d) is the average of (a) and (b). 

error (d) = [/(0) + /(0) + hf(0) + \h 2 f"(0) + ■ ■ ■ ] \ 

-[hf(0)+ 1 -h 2 f(0)+ 1 -hy"(0) + ...] 
\-l)h s f"(0) = ±h 3 f"(0) (11.18) 
Similarly, the error in (c) is 

error (c) = h[f(0) + \hf'(0) + \h 2 f"{Q) + • • • ] 

-[hf(0) + lh 2 f'(0)+ 1 -h 2 f"(0) + ...] 
« -^f"(0) (11.19) 
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The errors in the (c) and (d) formulas are both therefore the order of h 3 . 

Notice that just as the errors in formulas (a) and (b) canceled to highest order when you averaged 
them, the same thing can happen between formulas (c) and (d). Here however you need a weighted 
average, with twice as much of (c) as of (d). [1/12 — 2/24 = 0] 

1(d) + |(c) = [f(x ) + f(x + h)]^ + f{xo + h/2) h (11.20) 

This is known as Simpson's rule. 
Simpson's Rule 

Before applying this last result, I'll go back and derive it in a more systematic way, putting it into the 
form you'll see most often. 

Integrate Taylor's expansion over a symmetric domain to simplify the algebra: 

f h dx f{x) = 2/i/(0) + h 3 f" (0) + ^/""(O) + • • • (11.21) 

I'll try to approximate this by a three point formula af(—h) + /3f(0) + n ff(h) where a, f3, and 7, are 
unknown. Because of the symmetry of the problem, you can anticipate that a = 7, but let that go for 
now and it will come out of the algebra. 

af(-h)+pf(0) + lf(h) = 

a[f(0) - hf'(0) + ^ 2 /"(0) " \h s f'(0) + ^V""(0) + • • • ] 
+0/(0) 

+7[/(0) + hf'(0) + \tff> (0) + ±h 3 f"\0) + ^h*f""(0) + ■■■] 

You now determine the three constants by requiring that the two series for the same integral 
agree to as high an order as is possible for any f. 

2h = a + (5 + 7 
= —ah + 7/1 

\h 3 = \{a + 1 )h* 

and so, J* dx f{x) « | [f(-h) + 4/(0) + f(h)] . (11.22) 
The error term (the "truncation error") is 

^[fi-h) + AfiO) + fi-h)]-J h h dxfix) n±.±h*f"\0)-±h*f""(0) = ±h 5 f""i0) (11.23) 

Simpson's rule is exact up through cubics, because the fourth and higher derivatives vanish in 
that case. It's worth noting that there is also an elementary derivation of Simpson's rule: Given three 
points, there is a unique quadratic in x that passes through all of them. Take the three points to be 
( — h, f(—h)) , (0, /(0)) , and (h, fih)) , then integrate the resulting polynomial. Express your answer 
for the integral in terms of the values of / at the three points, and you get the above Simpson's rule. 
This has the drawback that it gives no estimate of the error. 



a = 7 = h/3, 



4/i/3 
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To apply Simpson's rule, it is necessary to divide the region of integration into an even number 
of pieces and apply the above formula to each pair. 



J a 



dx f(x) « ^ [f(x ) + 4/Orx) + f(x 2 )] + ± [f(x 2 ) + 4/(x 3 ) + /(x 4 )] + 



+ | [/(sjv- 2 ) + 4/(xjv_i) + /(xjv)] 
[/(xo) + 4/(:n) + 2/(x 2 ) + 4/(x 3 ) + • • • + 4/(^-0 + f(x N )] (11.24) 



Example: 



lo 1 + ^ 

Divide the interval to 1 into four pieces, then 



zdx = 4 tan 1 x 



7Y 



„ 1 + x 2C?:r ~ 12 



1 + 4 



1 



+ 2 



+ 



1 + (1/4) 2 1 + (1/2) 2 1 + (3/4)2 1 + 1 



3.1415686 



as compared to tt = 3.1415927. . .. 

When the function to be integrated is smooth, this gives very accurate results. 

Integration 

If the integrand is known at all points of the interval and not just at discrete locations as for tabulated 
or experimental data, there is more freedom that you can use to gain higher accuracy even if you use 
only a two point formula: 

f(x)dx^a[f((3) + f(-P)} 

-h 

I could try picking two arbitrary points, not symmetrically placed in the interval, but the previous 
experience with Simpson's rule indicates that the result will come out as indicated. (Though it's easy 
to check what happens if you pick two general points in the interval.) 

2/i/(0) + \h?f"{<S) + ^h*f»"(0) + ■■■ = a [2/(0) + /3 2 r (0) + ±pf""{0) + ■ ■ ■] 



To make this an equality through the low orders implies 



2h = 2a 
a = h 



with an error term 



-h 3 = a(3 2 

p = h/Vs 
i 



- ■ -h 5 f""(0) - -h 5 f""(0) = --^-h 5 f""(0), 
12 9 60 135 J w 



and 



(11.25) 



f f(x) dx^h[f (h/V3) + f (-h/V3)] [ - tU 5 /""(0)] (11-26) 
J — h 



With just two points, this expression yields an accuracy equal to the three point Simpson formula. 
Notice that the two points found in this way are roots of a certain quadratic 
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which is proportional to 



3 a; 2 - 1 = P2(x) 
2 2 y ' 



(11.27) 



the Legendre polynomial of second order. Recall section 4.11. 

This approach to integration, called Gaussian integration, or Gaussian quadrature, can be ex- 
tended to more points, as for example 



fix) dxnafi-P) +7/(0) + af(/3) 



The same expansion procedure leads to the result 



h 
9 



^[-!'\ll)+mo)+f(hJ 3 5 



wi 



(11.28) 



th an error proportional to h 7 f( 6 \0). The polynomial with roots 0, ±-y/3/5 is 



X 



x = P 3 (x) 



(11.29) 



the third order Legendre polynomial. 

For an integral / f(x) dx, let x = \{a + b)/2] + z[(b — a)/2\. Then for the domain — 1 < z < 1, 
x covers the whole integration interval. 

'' h b-a f 1 

f(x)dx= — - — / dzf{x) 



When you use an integration scheme such as Gauss's, it is in the form of a weighted sum over points. 
The weights and the points are defined by equations such as (11.26) or (11.28). 



J dzf(z) -> J2w k f{z k ) 



(11.30) 



or 



f(x) dx 



-i 
b — a 



£ w kf(%k), %k = [(a + b)/2] + z k [(b - a)/2) 

k 



Many other properties of Gaussian integration are discussed in the two books by C. Lanczos, 
"Linear Differential Operators," "Applied Analysis," both available in Dover reprints. The general 
expressions for the integration points as roots of Legendre polynomials and expressions for the coeffi- 
cients are there. The important technical distinction he points out between the Gaussian method and 
generalizations of Simpson's rule involving more points is in the divergences for large numbers of points. 
Gauss's method does not suffer from this defect. In practice, there is rarely any problem with using 
the ordinary Simpson rule as indicated above, though it will require more points than the more elegant 
Gauss's method. When problems do arise with either of these methods, they often occur because the 
function is ill-behaved, and the high derivatives are very large. In this case it can be more accurate 
to use a method with a lower order derivative for the truncation error. In an extreme case such a 
integrating a non-differentiable function, the apparently worst method, Eq. (11.15)(a), can be the best. 
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11.5 Differential Equations 

To solve the first order differential equation 

y' = f(x,y) y(x )=y , (11.31) 

the simplest algorithm is Euler's method. The initial conditions are y(xo) = yo, and y'(xo) = f(xo, yo), 
and a straight line extrapolation is 

y(x + h) =y + hf(x ,y ). (11.32) 

You can now iterate this procedure using this newly found value of y as a new starting condition to go 
from xq + h to Xq + 2h. 

Runge-Kutta 

Euler's method is not very accurate. For an improvement, change from a straight line extrapolation 
to a parabolic one. Take Xq = to keep the algebra down, and try a solution near in the form 
y{x) = a + f3x + 7a; 2 ; evaluate a, f3, and 7 so that the differential equation is satisfied near x = 0, 

y' = (3 + 2*yx = f(x, a + (3x + 7a; 2 ). 

Recall the Taylor series expansion for a function of two variables, section 2.5: 

f(x, y) = f(x , yo) + (x- x )D l f(x , y )+{y - y )D 2 f(x , y ) + - (x - xofDiDJixo, yo) 

+ \{y - yo) 2 D2D 2 f(x ,y ) + (x- x )(y - y )D 1 D 2 f(x , y ) + ■ ■ ■ (11.33) 

=> f3 + 2 1X = /(0, a) + xDJ(0, a) + (f3x + -fX 2 )D 2 f(0, «) + ■••. (11.34) 

The initial condition is at x = 0, y = yo, so a = yo- Equate coefficients of powers of x as high as is 
possible (here through rr 1 ). 

= /(0, a) 2 7 = D 1 f(0, a) + (3D 2 f(0, a). 

(If you set 7 = 0, this is Euler's method.) 

y(h) =y + hf(0, y ) + y [A/(0, yo) + f(0, y )D 2 f(0, y )] . (11.35) 

The next problem is to evaluate these derivatives. If you can easily do them analytically, you can 
choose to do that. Otherwise, since they appear in a term that is multiplied by ti 2 , it is enough to use 
the simplest approximation for the numerical derivative, 

Dif(0, Vo) = [m Vo) - /(0, yo)] /h (11.36) 

You cannot expect to use the same interval, h, for the y variable — it might not even have the same 
dimensions, 

D 2 f(0,yo)=[f(j,y + k)-fU,yo)]/k. (11.37) 

where j and k are the order of h. Note that because this term appears in an expression multiplied by 
h? , it doesn't matter what j is. You can choose it for convenience. Possible values for these are 

(1) J = k = hf(0,y o ) (3)j = h k = hf(0,y o ) 

(2) j = k = hf(h, y ) (4)j = h k = hf(h, y Q ). 
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The third of these choices for example gives 

h 2 

y = yo + hf(0,y ) + — 
= Vo + 3/(0,2/0) + ^f(h,y + hf(0,yo)) (11.38) 

This procedure, a second order Runge-Kutta method, is a moderately accurate method for advancing 
from one point to the next in the solution of a differential equation. It requires evaluating the function 
twice for each step of the iteration. 
Example: y' = 1 + y 2 y(0) = 0. Let h=0.1 



i[/(^o)-/(0, y o)] + /(0,,o) ^|^ 



X 


y(11.32) 


y(11.38) 


y(11.41) 


tan x 


0. 


0. 


0. 


0. 


0. 


0.1 


0.10 


0.10050 


0.10025 


0.10053 


0.2 


0.20100 


0.20304 


0.20252 


0.20271 


0.3 


0.30504 


0.30981 


0.30900 


0.30934 


0.4 


0.41435 


0.42341 


0.42224 


0.42279 


0.5 


0.53151 


0.54702 


0.54539 


0.54630 



(11.39) 



The fractional error at X = 0.5 with the second order method of equation (11.38) is 0.13% and with 
Euler it is —2.7% (first column). The method in Eq. (11.41) has error —0.17%. The fourth order 
version of Eq. (11.42) has an error —3.3 x 10~ 7 , more places than are shown in the table. 

The equation (11.38) is essentially a trapezoidal integration, like Eq. (11.15) (d). You would like 
to evaluate the function y at x = h by 

y(h)=y + / dxf(x,y(x)) (11.40) 
J 

A trapezoidal integration would be 

y(h)^yo + ^[f{o,y(o)) + f{h,y(h))] 

BUT, you don't know the y(h) that you need to evaluate the last term, so you estimate it by using the 
approximate result that you would have gotten from the simple Euler method, Eq. (11.32). The result 
is exactly equation (11.38). 

If you can interpret this Runge-Kutta method as just an application of trapezoidal integration, 
why not use one of the other integration methods, such as the midpoint rule, Eq. (11.15) (c)? This 
says that you would estimate the integral in Eq. (11.40) by estimating the value of / at the midpoint 
< x < h. 

y(h) ^y + hf{§, </(f )) ^y + y + |/(0, y )) (11.41) 

This is the same order of accuracy as the preceding method, and it takes the same number of function 
evaluations as that one (two). 

A commonly used version of this is the fourth order Runge-Kutta method: 



V = Vo + r [*1 + 2/c 2 + 2/c 3 + h] 



(11.42) 
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k 1 = hf(0,y o ) k 2 = hf{h/2,y + k 1 /2) 

h = hf (h/2, y + h/2) h = hf (h, y + k 3 ) ■ 

You can look up a fancier version of this called the Runge-Kutta-Fehlberg method. It's one of the 
better techniques around. 

Higher Order Equations 

How can you use either the Euler or the Runge-Kutta method to solve a second order differential 
equation? Answer: Turn it into a pair of first order equations. 

V" = f(x,y,y') — >y' = v, and v' = f(x, y, v) (11.43) 

The Euler method, Eq. (11.32) becomes 

y(x + h) =y(x ) + hv(x ), and v(x + h) = v (x ) + hf(x , y(x ), v{x )) 

The construction for Runge-Kutta is essentially the same. 
Adams Methods 

The Runge-Kutta algorithm has the advantage that it is self-starting; it requires only the initial condition 
to go on to the next step. It has the disadvantage that it is inefficient. In going from one step to the 
next, it ignores all the information available from any previous steps. The opposite approach leads to 
the Adams methods, though these are not as commonly used any more. I'm going to develop a little 
of the subject mostly to show that the methods that I've used so far can lead to disaster if you're not 
careful. 

Shift the origin to be the point at which you want the new value of y. Assume that you already 
know y at —h, —2h, . . . , —Nh. Because of the differential equation y' = f(x,y), you also know y' 
at these points. 

Assume 

N N 
3/(0) = £ a k y(-kh) + £ (3 k y'(-kh) (11.44) 
l l 

With 2N parameters, you can get this accurate to order h 21 ^^ 1 , 

y(-kh) = Y / (-kh) ny —P 
o n - 

Substitute this into the equation for y(0): 

y(o) = X> E(-^) n ^ + h f; h ±(-khr y( -^ 

k=l n=0 ' k=l n=0 

This should be an identity to as high an order as possible. The coefficient of h° gives 



N 

1 = £«fc 
k=l 



(11.45) 
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The next orders are 

o = X>fc(-wo + >*X>fc 

k k 



(11.46) 



iV = 1 is Euler's method again. 
N = 2 gives 



Oil + «2 = 1 «1 + 2«2 = (3l + 02 

an + 4a 2 = 2(/3i + 2/3 2 ) on + 8a 2 = 3(/?i + 4/3 2 ) 

The solution of these equations is 

ai = -4 a 2 = +5 ft = +4 /5 2 = +2 

j/(0) = -4j/(-/i) + 5y(-2h) + h[4y'(-h) + 2y'(-2h)} (11.47) 

To start this algorithm off, you need two pieces of information: the values of y at —h and at —2h. 
This is in contrast to Runge-Kutta, which needs only one point. 

Example: Solve y' = y y(0) = 1 (h = 0.1) 
I could use Runge-Kutta to start and then switch to Adams as soon as possible. For the purpose of 
this example, I'll just take the exact value of y at x = 0.1. 

e' 1 = 1.105170918 
y(.2) = -4y(.l) + 5y(0) + .1 [4/(.l, y(.l)) + 2/(0, 3/(0))] 
= -4y(.l) + 5y(0) + .4y(.l) + .2y(0) 
= -3.6y(.l) + 5.2y(0) 
= 1.221384695 

The exact value is e 2 = 1.221402758; the first error is in the underlined term. Continuing the calculation 
to higher values of x, 



X 


y 


.3 


1.3499038 


.4 


1.491547 


.5 


1.648931 


.6 


1.81988 


.7 


2.0228 


.8 


2.1812 


.9 


2.666 


1.0 


1.74 


1.1 


7.59 


1.2 


-18.26 


1.3 


105.22 



0. .5 1. 

Everything is going very smoothly for a while, though the error is creeping up. At around x = 1, 
the numerical solution goes into wild oscillation and is completely unstable. The reason for this is in 



11 — Numerical Analysis 



15 



the coefficients —4 and +5 of y(—h) and y(—2h). Small errors are magnified by these large factors. 
(The coefficients of y' are not any trouble because of the factor h in front.) 

Instability 

You can compute the growth of this error explicitly in this simple example. The equation (11.47) 
together with y' = y is 

y(0) = -3.6y(-h) + 5.2y(-2h), 

or in terms of an index notation 

y n = -3.6y n _i + 5.2|/ n _ 2 

This is a linear, constant coefficient, difference equation, and the method for solving it is essentially 
the same as for a linear differential equation — assume an exponential form y n = k n . 

k n = -3.6k n ^ + 5.2k n ~ 2 

k 2 + 3.6k -5.2 = 
k = 1.11 and -4.71 

Just as with the differential equation, the general solution is a linear combination of these two functions 
of n: 

y n = A(\.ll) n + B{-A.ll) n , 

where A and B are determined by two conditions, typically specifying y\ and yi. If B = 0, then y n 
is proportional to 1.11™ and it is the well behaved exponential solution that you expect. If, however, 
there is even a little bit of B present (perhaps because of roundoff errors), that term will eventually 
dominate and cause the large oscillations. If B is as small as 10~ 6 , then when n = 9 the unwanted 
term is greater than 1. 

When I worked out the coefficients in Eq. (11.47) the manipulations didn't look all that different 
from those leading to numerical derivatives or integrals, but the result was useless. This is a warning. 
You're in treacherous territory here; tread cautiously. 

Are Adams-type methods useless? No, but you have to modify the development in order to get 
a stable algorithm. The difficulty in assuming the form 

N N 

y(0) = £ a k y(-kh) + £ (3 k y'{-kh) 
i i 

is that the coefficients are too large. To cure this, you can give up some of the 2N degrees of 
freedom that the method started with, and pick the a priori to avoid instability. There are two 
common ways to do this, consistent with the constraint that must be kept on the a's, 

TV 

One way is to pick all the ot^ to equal l/JV. Another way is to pick a% = 1 and all the others = 0, 
and both of these methods are numerically stable. The book by Lanczos in the bibliography goes into 
these techniques, and there are tabulations of these and other methods in Abramowitz and Stegun. 

Backwards Iteration 

Before leaving the subject, there is one more kind of instability that you can encounter. If you try to 
solve y" = +y with y(0) = 1 and y'(0) = —1, the solution is e~ x . If you use any stable numerical 
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algorithm to solve this problem, it will soon deviate arbitrarily far from the desired one. The reason is 
that the general solution of this equation is y = Ae x + Be~ x . Any numerical method will, through 
rounding errors, generate a little bit of the undesired solution, e +x . Eventually, this must overwhelm 
the correct solution. No algorithm, no matter how stable, can get around this. 

There is a clever trick that sometimes works in cases like this: backwards iteration. Instead of 
going from zero up, start at some large value of x and iterate downward. In this direction it is the 
desired solution, e~ x , that is unstable, and the e +x is damped out. Pick an arbitrary value, say x = 10, 
and assign an arbitrary value to 2/(10), say 0. Next, pick an arbitrary value for ?/(10), say 1. Use these 
as initial conditions (terminal conditions?) and solve the differential equation moving left; necessarily 
the dominant term will be the unstable one, e~ x , and independent of the choice of initial conditions, 
it will be the solution. At the end it is only necessary to multiply all the terms by a scale factor to 
reduce the value at x = to the desired one; automatically, the value of y'(0) will be correct. What 
you are really doing by this method is to replace the initial value problem by a two point boundary value 
problem. You require that the function approach zero for large x. 

11.6 Fitting of Data 

If you have a set of data in the form of independent and dependent variables (i = 1, . . . ,N), 

and you have proposed a model that this data is to be represented by a linear combination of some set 
of functions, fn(x) 



M 



y = J2 a M x ^ 



(11.48) 



what values of will represent the observations in the best way? There are several answers to this 
question depending on the meaning of the word "best." The most commonly used one, largely because 
of its simplicity, is Gauss's method of least squares. 

Here there are N data and there are M functions that I will use to fit the data. You have to 
pick the functions for yourself. You can choose them because they are the implications of a theoretical 
calculation; you can choose them because they are simple; you can choose them because your daily 
horoscope suggested them. The sum of functions, X^m/a*' now depends on only the M parameters 
a^. The fs are fixed. The difference between this sum and the data points yi is what you want to be 
as small as possible. You can't use the differences themselves because they will as likely be negative as 
positive. The least squares method uses the sum of the squares of the differences between your sum of 
functions and the data. This criterion for best fit is that the sum 

2 

= Na 2 (11.49) 

be a minimum. The mean square deviation of the theory from the experiment is to be least. This 
quantity a 2 is called the variance. 

Some observations to make here: iV > M, for otherwise there are more free parameters than 
data to fit them, and almost any theory with enough parameters can be forced to fit any data. Also, 
the functions must be linearly independent; if not, then you can throw away some and not alter the 



N 



M 
H=l 
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result — the solution is not unique. A further point: there is no requirement that all of the Xj are 
different; you may have repeated the measurements at some points. 

Minimizing this is now a problem in ordinary calculus with M variables. 



d 



N 

£ 

i=i 



M 
//=! 



UiXi) = 



rearrange: 



a. 



These linear equations are easily expressed in terms of matrices. 

Ca = b, 



where 



N 



i=l 



(11.50) 



(11.51) 



a is the column matrix with components and b has components YliVifv{ x i)- 
The solution for a is 

a = C~ 1 b. (11.52) 

If C turned out singular, so this inversion is impossible, the functions were not independent. 
Example: Fit to a straight line 



Then Ca = b 



is 



hix) = 1 



N ^xA fax 
E%i Erf J V«2 



hix) = x 



E Vi x i 



The inverse is 



a i 



E*i -E^A/Ei/: 



and the best fit line is 

y = ax + a^s 



^ iV 



(11.53) 



11.7 Euclidean Fit 

In fitting data to a combination of functions, the least squares method used Eq. (11.49) as a measure 
of how far the proposed function is from the data. If you're fitting data to a straight line (or plane if 
you have more variables) there's another way to picture the distance. Instead of measuring the distance 
from a point to the curve vertically using only y, measure it as the perpendicular distance to the line. 
Why should this be any better? It's not, but it does have different uses, and a primary one is data 
compression. 
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y 




Do this in two dimensions, fitting the given data to a straight line, and to describe the line I'll 
use vector notation, where the line is u + av and the parameter a varies over the reals. First I need 
to answer the simple question: what is the distance from a point to a line? The perpendicular distance 
from w to this line requires that 

d 2 



[w — u — av ) 

be a minimum. Differentiate this with respect to a and you have 



(w — u — av) ■ ( — v) = implying 
For this value of a what is d 2 7 

d 2 



av 



(w — u) -v 



(w - u) 2 + a 2 v 2 



2av- (w — u) 



{w-uf - ^[(w-u)-v] 2 



(11.54) 



Is this plausible? (1) It's independent of the size of v, depending on its direction only. (2) It depends 
on only the difference vector between w and u, not on any other aspect of the vectors. (3) If I add any 
multiple of v to u, the result is unchanged. See problem 11.37. Also, can you find an easier way to get 
the result? Perhaps one that simply requires some geometric insight? 

The data that I'm trying to fit will be described by a set of vectors Wi, and the sum of the 
distances squared to the line is 



D 2 



N 



U 



2 N 1 



Now to minimize this among all u and v I'll first take advantage of some of the observations from the 
preceding paragraph. Because the magnitude of v does not matter, I'll make it a unit vector. 



D2 = E 



(11.55) 



Now to figure out u: Note that I expect the best fit line to go somewhere through the middle of the 
set of data points, so move the origin to the "center of mass" of the points. 



w n 



and let 



w; 



Wi - w n 



and 



then the sum = ar| d 

D 2 = Y,w' 2 + Nu 12 ~Y,{wl-v) 2 - N(u' -v)'- 



u 



u — w n 



(11.56) 



This depends on four variables, u' x , u' y , v x and v y . If I have to do derivatives with respect to all of them, 
so be it, but maybe some geometric insight will simplify the calculation. I can still add any multiple of v 
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to u without changing this expression. That means that for a given v the derivative of D 2 as u' changes 
in that particular direction is zero. It's only as u' changes perpendicular to the direction of v that D 2 
changes. The second and fourth term involve u' 2 — (u' ■ v ) 2 = u' 2 (l — cos 2 6) = u' 2 sin 2 6, where this 
angle 6 is the angle between u' and v. This is the perpendicular distance to the line (squared). Call it 



u'j_ = u' sm.6. 



D 2 = Y, w? ~ J2(Wi ■ v? + Nu' 2 - N{u' -v) 2 = Y^ w i ~ ■ v? + Nu'l 

The minimum of this obviously occurs for u'j_ = 0. Also, because the component of u' along the 
direction of v is arbitrary, I may as well take it to be zero. That makes u' = 0. Remember now that 
this is for the shifted w ' data. For the original Wi data, u is shifted to u = w n 



'mean ■ 



I'm not done. What is the direction of v? That is, I have to find the minimum of D 2 subject to 
the constraint that \v\ = 1. Use Lagrange multipliers (section 8.12). 

Minimize D 2 = V] w' 2 - ^(w- ■ v ) 2 subject to (f> = vl + v 2 -l = 
The independent variables are v x and v y , and the problem becomes 

V (D 2 + \<j>) = 0, with = 

Differentiate with respect to the independent variables and you have linear equations for v x and Vy, 

5 2 -}22(w' xi v x +w' yi v y )w xi + \2v x = 

"ArLKi^ + W +\2v x = or 

(11.58) 

Correlation, Principal Components 

The correlation matrix of this data is 



(C) 



N \ Yw' -w' ■ y w 12 - 



The equations (11.58) are 



C XX C X y\ /V X \ = y fv X \ (1L5g) 

^yx <^yy J \ v y J \ V V J 

where A' = X/N. This is a traditional eigenvector equation, and there is a non-zero solution only if the 
determinant of the coefficients equals zero. Which eigenvalue to pick? There are two of them, and one 
will give the best fit while the other gives the worst fit. Just because the first derivative is zero doesn't 
mean you have a minimum of D 2 ; it could be a maximum or a saddle. Here the answer is that you 
pick the largest eigenvalue. You can see why this is plausible by looking at the special case for which 
all the data lie along the x-axis, then C xx > and all the other components of the matrix = 0. The 
eigenvalues are C xx and zero, and the corresponding eigenvectors are x and y respectively. Clearly the 
best fit corresponds to the former, and the best fit line is the x-axis. The general form of the best fit 
line is (now using the original coordinate system for the data) 

1 - . 

av + AT 2^ Wi = av + Wmean 
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and this v is the eigenvector having the largest eigenvalue. More generally, look at Eq. (11.57) and you 
see that the lone negative term is biggest if the ws are in the same direction (or opposite) as v. 

y vA 




This establishes the best fit to the line in the Euclidean sense. What good is it? It leads into 
the subject of Principal Component Analysis and of Data Reduction. The basic idea of this scheme is 
that if this fit is a good one, and the original points lie fairly close to the line that I've found, I can 
replace the original data with the points on this line. The nine points in this figure require 9 x 2 = 18 
coordinates to describe their positions. The nine points that approximate the data, but that lie on the 
line and are closest to the original points require 9x1 = 9 coordinates along this line. Of course you 
have some overhead in the data storage because you need to know the line. That takes three more 
data (u and the angle of v), so the total data storage is 12 numbers. See problem 11.38 

This doesn't look like much of a saving, but if you have 10 6 points you go from 2 000 000 
numbers to 1000 003 numbers, and that starts to be significant. Remember too that this is a two 
dimensional problem, with two numbers for each point. With more coordinates you will sometimes 
achieve far greater savings. You can easily establish the equation to solve for the values of a for each 
point, problem 11.38. The result is 

ttj = (wi — u) • v 



11.8 Differentiating noisy data 

Differentiation involves dividing a small number by another small number. Any errors in the numerator 
will be magnified by this process, and if you have to differentiate experimental data this will always 
present a difficulty. If it is data from the output of a Monte Carlo calculation the same problem will 
arise. 

Here is a method for differentiation that minimizes the sensitivity of the result to the errors in 
the input. Assume equally spaced data where each value of the dependent variable f(x) is a random 
variable with mean (/(x)) and variance a 2 . Follow the procedure for differentiating smooth data and 
expand in a power series. Let h = 2k and obtain the derivative between data points. 

f(k) = /(o) + kf(o) + \k 2 f{Q) + + ■■■ 

f(k)-f(-k) = 2kf(0)+ 1 -k^'"(0) + ... 
f(3k) - f(-3k) = 6*/'(0) + y k 3 f"'(0) + ■■■ 
I'll seek a formula of the form 



/' (0) = a [/(A;) - /(-*)] +(3[f(3k) - f(-3k)] (11.60) 

I am assuming that the variance of / at each point is the same, a 2 , and that the fluctuations in / at 
different points are uncorrelated. The last statement is, for random variables f± and / 2 , 



<(/i ~ </i» (/a - </ 2 »> = which expands to </i/ 2 ) = </i)</ 2 > (11.61) 
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Insert the preceding series expansions into Eq. (11.60) and match the coefficients of /'(0). This 
gives an equation for a and (3: 

2ka + 6k/3 = l (11.62) 

One way to obtain another equation for a and (3 is to require that the k 3 f"'(0) term vanish; this leads 
back to the old formulas for differentiation, Eq. (11.11). Instead, require that the variance of /'(0) be 
a minimum. 

<(f (0) - <f (0)» 2 > = ([<*(/(*) - (/(*)» +a(f(-k) - (f(-k))) + • • -] 2 > 

= 2a 2 a 2 + 2o 2 f3 2 (11.63) 

This comes from the fact that the correlation between say f(k) and f(—3k) vanishes, and that all the 
individual variances are a 2 . That is, 



(/(*)-</(*)))(/(-*)-</(-*)))) = 



along with all the other cross terms. Problem: minimize 2a 2 (a 2 + /3 2 ) subject to the constraint 
2ka + 6k/3 = 1. It's hardly necessary to resort to Lagrange multipliers for this problem. 
Eliminate a: 



dj3 



(3 = 3/20A;, a = l/20fc 



fm M -3/(-ft)-/(»)+/W + 3/(2A) | (u M) 

and the variance is 2a 2 (a 2 + f3 2 ) = a 2 /5h 2 . In contrast, the formula for the variance in the standard 
four point differentiation formula Eq. (11.10), where the truncation error is least, is 65cr 2 /72h 2 , which 
is 4.5 times larger. 

When the data is noisy, and most data is, this expression will give much better results for this 
derivative. Can you do even better? Of course. You can for example go to higher order and both 
decrease the truncation error and minimize the statistical error. See problem 11.22. 

11.9 Partial Differential Equations 

I'll illustrate the ideas involved here and the difficulties that occur in even the simplest example of a 
PDE, a first order constant coefficient equation in one space dimension 

du/dt + cdu/dx = u t + cu x = 0, (11.65) 

where the subscript denotes differentiation with respect to the respective variables. This is a very simple 
sort of wave equation. Given the initial condition that at t = 0, u(0,x) = f(x), you can easily check 
that the solution is 

u(t,x) = f{x - ct) (11.66) 

The simplest scheme to carry data forward in time from the initial values is a generalization of 
Euler's method for ordinary differential equations 

u(t + At, x) = u(t, x) + Ut(t, x)At 
= u(t, x) — u x (t, x)cAt 
cAt 

= u(t, x) - — — [u(t, x + Ax) - u(t, x - Ax)] , (11.67) 
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Here, to evaluate the derivative, I used the two point differentiation formula. 

• ••••• 

In this equation, the value of u at point (At, 4Ax) depends on the values at (0, 3 Ax), (0, 4Ax), 
and (0,5 Ax). This diagram shows the scheme as a picture, with the horizontal axis being x and the 
vertical axis t. You march the values of u at the grid points forward in time (or backward) by a set of 
simple equations. 

The difficulties in this method are the usual errors, and more importantly, the instabilities that 
can occur. The errors due to the approximations involved can be classified in this case by how they 
manifest themselves on wavelike solutions. They can lead to dispersion or dissipation. 

Analyze the dispersion first. Take as initial data u(t,x) = Acoskx (or if you prefer, e ikx ). The 
exact solution will be Acos(kx — cut) where u = ck. Now analyze the effect of the numerical scheme. 
If Ax is very small, using the discrete values of Ax in the iteration give an approximate equation 



2Ax 



[u(t, x + Ax) — u(t, x — Ax)] 



A power series expansion in Ax gives, for the first two non-vanishing terms 



u x + \(Ax) 2 u 2 

D 



(11.68) 



So, though I started off solving one equation, the numerical method more nearly represents quite a 
different equation. Try a solution of the form Acos(kx — ut) in this equation and you get 



u = c 



k--(Ax) 2 k 3 



(11.69) 



and you have dispersion of the wave. The velocity of the wave, u/k, depends on k and so it depends 
on its wavelength or frequency. 

The problem of instabilities is more conveniently analyzed by the use of an initial condition 



u(0,x) = e ikx , then Eq. (11.67) is 



u(At, x) = e 



ikx 



cAt 



Ak(x+Ax) _ ik(x—Ax) 



Akx 



2Ax 
icAt 



1 



Ax 



sin 



kAx 



(11.70) 



The n-fold iteration of this, therefore involves just the n th power of the bracketed expression; that's 
why the exponential form is easier to use in this case. If A;Ax is small, the first term in the expansion 
of the sine says that this is approximately 

e ikx [l-ikcAt] n , 
and with small At and n = t/ At a large number, this is 



Akx 



ikct 



Ak(x—ct) 



11 
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1 C 2 (At) 2 . 9 i a 

1 + ' sm 2 kAx 
(Ax) 2 



1/2 



> 1 



(11.71) 



so the magnitude of the solution grows exponentially. This instability can be pictured as a kind of 
negative dissipation. This growth is reduced by requiring kcAt -C 1. 

Given a finite fixed time interval, is it possible to get there with arbitrary accuracy by making At 
small enough? With n steps = t/ At, r n is 



1 c 2 (At) 
1 + jAxY" ] 

[1 + a] 1 /" 

cHAt 



2 kAx 



t/2At 



[1 + a] 



exp 



2(Ax) ; 



sin 2 kAx 



so by shrinking At sufficiently, this is arbitrarily close to one. 

There are several methods to avoid some of these difficulties. One is the Lax-Friedrichs method: 

u(t + At, x) = i [u(t, x + Ax) + u(t, x - Ax)] - |^ [u(t, x + Ax) - u(t, x - Ax)] (11.72) 

By appropriate choice of At and Ax, this will have r < 1, causing a dissipation of the wave. Another 
scheme is the Lax- Wend roff method. 



u(t + At, x) = u(t, x) 



cAt 
2Ax 
c 2 (At) 2 



2(Ax) 

This keeps one more term in the power series expansion. 



[u(t, x + Ax) — u(t, x — Ax)] 
., [u(t, x + Ax) — 2u(t, x) + u(t, x — Ax)] 



(11.73) 



Exercises 

1 Use the four-point interpolation formula Eq. (11.3) to estimate e 3 / 4 from the values of e x at 0, 1/2, 
1, 3/2. From the known value of the number, compute the relative error. 

2 Find a root of the equation cosx = x. Start with a graph of course. 

3 Find the values of a and of x for which e x = ax has a single root. 

4 Find the roots of e x = ax when a is twice the value found in the preceding exercise, (and where is 
your graph?) 

5 Use (a) midpoint formula and (b) Simpson's rule, with two intervals in each case, to evaluate 

4 /(J dx 1/(1 + x 2 ). 

6 Use Newton's method to solve the equation sinx = 0, starting with the initial guess Xq = 3. 
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Problems 

11.1 Show that a two point extrapolation formula is 

f(0)^2f(-h)-f(-2h) + h 2 f"(0). 

11.2 Show that a three point extrapolation formula is 

/(0) « 3f(-h) - 3f(-2h) + f(-3h) + h 3 f"'(0). 

11.3 Solve x 2 — a = by Newton's method, showing graphically that in this case, no matter what the 
initial guess is (positive or negative), the sequence will always converge. Draw graphs. Find (This 
is the basis for the library square root algorithm on some computers.) 

11.4 Find all real roots of e~ x = sinx to ±1(T 4 . Ans: 0.588533, vr - 0.045166, 2vr + 0.00187 . . . 

11.5 The first root r\ of e~ ax = s'mx is a function of the variable a > 0. Find dr\jda at a = 1 by 
two means, (a) First find r\ for some values of a near 1 and use a four-point differentiation formula, 
(b) Second, use analytical techniques on the equation to solve for dr\/da and evaluate the derivative 
in terms of the known value of the root from the previous problem. 

11.6 Evaluate erf(l) = ^= j Q l dte~ t2 Ans: 0.842736 (more exact: 0.842700792949715) 

11.7 The principal value of an integral is (a < Xq < b) 



p f b Jix)_ dx = Hm 

J a X-X Q 



a X Xq Jxo+e x x 



(a) Show that an equal spaced integration scheme to evaluate such an integral is (using points 0, ±h) 

"+h 
-h 

(b) Also, an integration scheme of the Gaussian type is 



' J +h h ^ dx = f(h) - n-h) - 2 -h*r(o). 



V3[f(h/Vs) - f(-h/Vs)] + §- 5 f v (o)- 



11.8 Devise a two point Gaussian integration with errors for the class of integrals 

/■+00 2 

/ dxe x f(x). 

Find what standard polynomial has roots at the points where / is to be evaluated. 

Ans:|V5F[/(-l/v^) + /(l/v^)] 

11.9 Same as the preceding problem, but make it a three point method. 
Ans: - iV§) + 1/(0) + |/( + |V5)] 
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11.10 Find one and two point Gauss methods for 

/•oo 

/ dxe~ x f(x). 
Jo 

(a) Solve the one point method completely. 

(b) For the two point case, the two points are roots of the equation 1 — 2x + \x 2 = 0. Use that as 
given to find the weights. Look up Laguerre polynomials. 

11.11 In numerical differentiation it is possible to choose the interval too small. Every computation is 
done to a finite precision, (a) Do the simplest numerical differentiation of some specific function and 
take smaller and smaller intervals. What happens when the interval gets very small? (b) To analyze 
the reason for this behavior, assume that every number in the two point differentiation formula is kept 
to a fixed number of significant figures (perhaps 7 or 8). How does the error vary with the interval? 
What interval gives the most accurate answer? Compare this theoretical answer with the experimental 
value found in the first part of the problem. 

11.12 Just as in the preceding problem, the same phenomenon caused by roundoff errors occurs in 
integration. For any of the integration schemes discussed here, analyze the dependence on the number 
of significant figures kept and determine the most accurate interval. (Surprise?) 

11.13 Compute the solution of y' = 1 + y 2 and check the numbers in the table where that example 
was given, Eq. (11.39). 

11.14 If in the least square fit to a linear combination of functions, the result is constrained to pass 
through one point, so that ^ a^f^Xo) = K is a requirement on the a's, show that the result becomes 

a = C- 1 [b + \f ], 

where fo is the vector /^(xo) and A satisfies 

A(/o,C- 1 /o)=^-</o,C- 1 6>. 



11.15 Find the variances in the formulas Eq. (11.8) and (11.10) for /', assuming noisy data. 
Ans: a 2 /2h 2 , 65a 2 /72h 2 

11.16 Derive Eqs. (11.61), (11.62), and (11.63). 

11.17 The Van der Pol equation arises in (among other places) nonlinear circuits and leads to self- 
exciting oscillations as in multi-vibrators 

d x , o . dx 

Take e = .3 and solve subject to any non-zero initial conditions. Solve over many periods to demonstrate 
the development of the oscillations. 

11.18 Find a formula for the numerical third derivative. Cf. Eq. (2.18) 

11.19 The equation resulting from the secant method, Eq. (11.7), can be simplified by placing every- 
thing over a common denominator, (f(x2) — f(x±)) . Explain why this is a bad thing to do, how it can 
lead to inaccuracies. 
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11.20 Rederive the first Gauss integration formula Eq. (11.25) without assuming the symmetry of the 
result 

r+h 

/ f(x)dx^af{(3) +1 f{5). 
J-h 

11.21 Derive the coefficients for the stable two-point Adams method. 

11.22 By putting in one more parameter in the differentiation algorithm for noisy data, it is possible 
both to minimize the variance in /' and to eliminate the error terms in h 2 f". Find such a 6-point 
formula for the derivatives halfway between data points OR one for the derivatives at the data points 
(with errors and variance). 

Ans: /'(0) = [58 (/(A) - f{-h)) + 67(/(2/i) - f(-2h)) - 22{f(3h) - f (-3hj)] / (252h) 

11.23 In the same spirit as the method for differentiating noisy data, how do you interpolate noisy 
data? That is, use some extra points to stabilize the interpolation against random variations in the 
data. To be specific, do a midpoint interpolation for equally spaced points. Compare the variance here 
to that in Eq. (11.3). Ans: f(0) « [f(-3k) + f(-k) + f (k) + /(3jfe)]/4, a 2 is 4.8 times smaller 

11.24 Find the dispersion resulting from the use of a four point formula for u x in the numerical solution 

of the PDE u t + cu x = 0. 

11.25 Find the exact dispersion resulting from the equation 

ut = —c[u(t, x + Ax) — u(t, x — Ax)]/ 2 Ax. 
That is, don't do the series expansion on Ax. 

11.26 Compute the dispersion and the dissipation in the Lax-Friedrichs and in the Lax- Wend roff meth- 
ods. 

11.27 In the simple iteration method of Eq. (11.71), if the grid points are denoted x = mAx, t = nAt, 
where n and m are integers (— oo < n,m < +oo), the result is a linear, constant-coefficient, partial 
difference equation. Solve subject to the initial condition 

u(0,m) = e ikmAx . 



11.28 Lobatto integration is like Gaussian integration, except that you require the end-points of the 
interval to be included in the sum. The interior points are left free. Three point Lobatto is the same 
as Simpson; find the four point Lobatto formula. The points found are roots of -P^-i- 

11.29 From the equation y' = f(x,y), one derives y" = f x + ffy. Derive a two point Adams type 
formula using the first and second derivatives, with error of order h 5 as for the standard four-point 
expression. This is useful when the analytic derivatives are easy. The form is 

y(0) = y(-h) + (3 iy '(-h) + (3 2 y'(-2h) + lx y"{-h) + j2y"{-2h) 

Ans: Pi = -h/2, (3 2 = 3h/2, 71 = 17/i 2 /l2, 72 = 7h 2 /l2 

11.30 Using the same idea as in the previous problem, find a differential equation solver in the spirit 
of the original Euler method, (11.32), but doing a parabolic extrapolation instead of a linear one. That 
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is, start from (xo,yo) and fit the initial data to y = a + (3{x — Xq) + 7(2; — Xo) 2 in order to take a 
step. Ans: y(h) = y + hf(0,y o ) + (/i 2 /2) [f x (0,y ) + / y (0,y o )/(0,y o )] 

11.31 Show that the root finding algorithm of Eq. (11.7) is valid for analytic functions of a complex 
variable with complex roots. 

11.32 In the Runge-Kutta method, pick one of the other choices to estimate /^/(O, yo) in Eq. (11.37). 
How many function evaluations will it require at each step? 

11.33 Sometimes you want an integral where the data is known outside the domain of integration. 
Find an integration scheme for f(x)dx in terms of f(h), /(0), and f(—h). Ans: [—/(—/&) + 
8/(0) + 5f(h))h/l2, error oc h 4 

11.34 When you must subtract two quantities that are almost the same size, you can find yourself 
trying to carry ridiculously many significant figures in intermediate steps. If a and b are very close and 
you want to evaluate \fa — Vb, devise an algorithm that does not necessitate carrying square roots out 
to many more places than you want in the final answer. Write a = b + e. 

Ans: e/2-A error: e 2 /8b 3 / 2 

11.35 Repeat the preceding problem but in a more symmetric fashion. Write a = x + e and b = x — e. 
Compare the sizes of the truncation errors. Ans: e/y/x, — e 3 /8x 5 / 2 

11.36 The value of tt was found in the notes by integrating 4/ (1+x 2 ) from zero to one using Simpson's 
rule and five points. Do the same calculation using Gaussian integration and two points. Ans: 3.14754 
(Three points give 3.14107) 

11.37 (a) Derive Eq. (11.54). 

(b) Explain why the plausibility arguments that follow it actually say something. 

11.38 After you've done the Euclidean fit of data to a straight line and you want to do the data 
reduction described after Eq. (11.59), you have to find the coordinate along the line of the best fit to 
each point. This is essentially the problem: Given the line [u and v) and a point (w), the new reduced 
coordinate is the a in u + av so that this point is closest to w. What is it? You can do this the hard 
way, with a lot of calculus and algebra, or you can draw a picture and write the answer down. 

11.39 Data is given as (X{,yi) = {(1,1), (2,2), (3,2)}. Compute the Euclidean best fit line. Also 
find the coordinates, «j, along this line and representing the reduced data set. 

Ans: u = (2,5/3) v = (0.88167, 0.47186) a x = -1.1962 a 2 = 0.1573 a 3 = 1.0390 
The approximate points are (0.945, 1.102), (2.139, 1.741), (2.916, 2.157) 

[It may not warrant this many significant figures, but it should make it easier to check your work.] 

11.40 In the paragraph immediately following Eq. (11.23) there's mention of an alternate way to derive 
Simpson's rule. Carry it out, though you already know the answer. 

11.41 (a) Derive the formula for the second derivative in terms of function values at three equally 
spaced points, (b) Use five points to get the second derivative, but using the extra data to minimize 
the sensitivity to noise in the data. Ans: (a) [f(-h) - 2/(0) + f{h)]/h 2 (b) [ - 2/(0) - (f(h) + 
f{-h)) + 2(/(2/i) + f(-2h))]/7h 2 . The ratio of the standard deviation for (b) is smaller than for (a) 
by a factor 21 1 / 2 = 4.6 



